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ABSTRACT 

In this paper we present a high-order Lagrangian-decoupling method for the unsteady 
convection-diffusion and incompressible Navier-Stokes equations. The method is based upon: 
Lagrangian variational forms that reduce the convection-diffusion equation to a symmetric 
initial value problem; implicit high-order backward-differentation finite- difference schemes 
for integration along characteristics; finite element or spectral element spatial discretiza- 
tions; mesh-invariance procedures and high-order explicit time-stepping schemes for deducing 
function values at convected space-time points. The method improves upon previous finite 
element characteristic methods through the systematic and efficient extension to high order 
accuracy, and the introduction of a simple structure-preserving characteristic-foot calcula- 
tion procedure which is readily implemented on modern architectures. The new method is 
significantly more efficient than explicit-convection schemes for the Navier-Stokes equations 
due to the decoupling of the convection and Stokes operators and the attendant jncrease^n 
temporal stability. Numerous numerical examples are given for the convection- diffusion and 
Navier-Stokes equations for the particular case of a spectral element spatial discretization. 
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1. Introduction 


A large class of important fluid flows is described by the incompressible Navier-Stokes 
equations, 


+ U 1 U m,9 = ~P,m + f m in n 

u 9lt = 0 in f2 , 


(la) 

(16) 


where x - (x m ) m =i d is the space variable in 12 of JR d , t is time, f_(x, t ) is the prescribed 

force, u(x,t) = (u m (x,t)) m=1 d is the velocity, p(x,t) is the pressure normalized by a 

(constant) density, and 1 / is the kinematic viscosity. We adopt Cartesian tensor indicial 
notation, with — 8s “ , _ — “ , and summation over repeated subscripts, u t y = 

^»=1 lit' Ec l uatlon ( la ) represents conservation of momentum for a Newtonian fluid, 
while equation (lb) represents conservation of mass for an essentially incompressible flow. 

The difficulties associated with numerical solution of the incompressible Navier-Stokes 
equations arise from two distinct sources: the pressure-divergence terms, and the con- 
vection diffusion contributions. The pressure-divergence terms require first, for optimal 


convergence, that proper discrete function spaces be chosen for the velocity and pressure 
[1], and second, for rapid solution, that algorithms be developed that effectively decouple 
the pressure and velocity [e.g., 2-4], The origin of the numerical difficulties associated 
with the convection and convection-diffusion terms are equally well documented, but less 
resolved: convection- diffusion introduces thin boundary layers for small u; the imbalance 
of convection and diffusion leads to weak spatial “stability” for small t/, in that the ratio 
of the continuity and coercivity constants becomes large [5]; timescales for equilibration 
become large for small t/; the convection term destroys the isotropy and symmetry of the 
Stokes operator for all values of u. 


To motivate our new class of schemes we briefly summarize standard numerical ap- 
proaches to the time-dependent Navier-Stokes equations. In particular, we note that the 
unsteady Navier-Stokes equations are often treated in a semi-implicit fashion, in which 
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the symmetric Stokes operator is handled implicitly and the convection term is treated 
explicitly, 

"m +> -J C + jr([u, Um ,,]",[ U ,u„,.]"- 1 ,...) = 

A t 

+ in ft (2a) 

U n+l = 0 in ft . ( 2b ) 

Here u" (x) = u(z, nAt), where At is the time step, and T and Q are appropriate explicit 
and implicit time-marching schemes [6], respectively. Explicit treatment of the convective 
term is motivated by the fact that solution of the (nonlinear) implicit convection equa- 
tions is difficult; no efficient (faster- than- time-like) robust iterative methods are available 
for general unsymmetric anisotropic operators, and direct methods are typically memory- 
intensive, costly, and serial, in particular in higher space dimensions. These arguments 
for explicit convection, although by no means universally accepted (e.g., (7)), are most 
easily defended for high-order spatial discretizations, in which asymmetry and bandwidth 
problems are further aggravated by longer-range coupling. 

The semi-implicit method of equation (2) is intended to represent a compromise of 
convenience and efficiency between the symmetric Stokes operator and the unsymmetric 
convection term; unfortunately, it is not an optimal compromise, as we now describe. To 
begin, we assume that the semi-discretization (2) is further discretized in space, with the 
spatial discretization characterized by an effective mesh spacing h. We next (plausibly) 
assume that the critical time step set by the explicit convection treatment, 

A t„ ~ 0(h/\u\), 

is smaller than the time step required for accuracy, At ace , and that stability thus determines 
the timestep. The work required to integrate (2) to a final time T, W* ‘- (semi-implicit), 

is then 

W*- = (T/At cr )WAi .( 3a ) 
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where Wa« is the total work per timestep (e.g., in .clock cycles). The total work per 
timestep Wa» is further broken up as: 


W At = + W£ t *, . 


( 36 ) 


Here the convective work per timestep is given by Wa*,* = aVV££, where W£ t a, e is the 
work xequired to evaluate the convective terms, and a is the (order unity) number of 
evaluations required for T , and the Stokes work per timestep, Wa* | 0 is equal to the work 
required to invert the Stokes system, 

Finally, we note that for nontrivial problems it is typically the case that W^ 1 , 0 ' << 
as ^ ie solution for the pressure in incompressible simulations requires at least an 
elliptic solve [8-10], and, more properly, nested elliptic solves [2-4] at each timestep. For 
instance, for typical spectral-element-discretization Navier-Stokes calculations, Wg t al is 
an order of magnitude less than W&*, for splitting Poisson schemes [11-13], and several 
orders-of-magnitude less than Wa", for consistent Uzawa methods [4,14,15], (Note that 
explicit treatment of the viscous terms does not significantly alter these estimates, as 
the majority of the work derives directly from the elliptic pressure term.) These 

complexity estimates are based on relatively efficient conjugate-gradient [14] and multigrid 
[16,17] elliptic solvers [4,14,15], and although 1V^ ( will certainly be reduced by further 

improvements in elliptic technology, it is clear that will always remain significantly 

larger than 

It is now readily seen in what sense semi-implicit methods are not entirely satisfactory. 
The usual motivation for turning to an explicit rather than implicit formulation is that 
the increase in number of timesteps (T j At) is more than counteracted by an associated 
decrease in work per timestep, Wa«. However, it is clear from (3) that for the Navier- 
Stokes equations an increase in T / At from T/ At ace to T / At er is not compensated for by 
a commensurate or greater decrease in Wa», as the majority of the work per timestep is 
due to the elliptic solves, Wa ( ,, >> WAt, c - A better compromise is required, in which 
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T/A( « T/A(„, and Wa,„ and Wa,., are of the same order: such a “balanced" scheme 
would clearly result iu improvements in computational efficiency, in particular for high- 
order spatial discretizations. These comments are especially relevant when (2) is viewed a 

method for the solution of the steady-state equations. 

We present in this paper a new liigh-order Lagrangian-decoupling method for the 
convection-diffusion and Navier-Stokes equations based on lugh-order finite-differences in 
“time" and variational (finite or spectral element) discretizations in space . The scheme 
preserves symmetry as regards the diffusion and Stokes operators, thereby allowing for fast 
and robust iterative solution, yet provides for much better temporal stability fo, the full 
Navier-Stokes equations than explicit-convection alternatives such as (2); in essence, the 
new method exchanges a significant increase in At above At„ for an easily accomodated 
increase in the convective work that brings Wa„ into hulnnc. with Wa,,- The direct an- 
cestors of the method are the Lagrangian and arbitrary-Lagrangian-Eulerian finite element 
techniques (18,19], however the technique can also be thought of as a rigorous subcycling 
approach (20), or as a hi igl, -order non-dissipative upwinding scheme [21,22], The essential 
differences between our method and its closest predecessor, the characteristic scheme of 
Pironneau [18], are the systematic and efficient extension to high-order accuracy, and a 
simple structure-preserving Lagrangian remesliing scheme (characteristic-foot calculation) 

which is readily implemented on modern architectures. 

The structure of the paper is as follows. In Section 2 we introduce the strong 
Lagrangian variational forms of the convection-diffusion equation. In Section 3 we describe 
the nodal function spaces for finite [23] or spectral [14,24] element spatial discretization; we 
present the new high-order fi„ite-differe.,ce-in-“time” Lagrangian-decoupling method; and 
we summarize the stability, accuracy, and computational complexity of the new technique. 
In Section 4 we briefly describe an extension of the scheme to the fuU unsteady nonlinear 
Navier-Stokes equations. Lastly, in Section 5, numerous numerical examples are given of 
the new scheme for the particular case of a spectral element spatial discretization. 
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2. Convection-Diffusion Equation 

2.1 Strong Form 

In tins paper we shall adopt the common practice of using the convection-diffusion 
equation as a simple model for the Navier-Stokes equations (1). As our scheme is partic- 
ularly tied to the structure of the incompressible Navier-Stokes equations (the inequality 
WAt.c << Wai,. is only occasionally true for convection-diffusion alone), it is imperative 
that we return to the parent equation to demonstrate that our method is, indeed, applica- 
ble to the full Navier-Stokes system. This extension is described and illustrated in Sections 
4 and 5, respectively. 

Our simple convection-diffusion equation is given by 

<t>« +u t <f>, t = i/0, „ + g vte(o,r), Vxefi , (4) 

where 4>{x,t) is a passive scalar, u is the prescribed convection velocity, x is space, t is 
time, t/ is diffusivity, g(x,t) is a prescribed source, and D is the domain in JR d . In addition 
to equation (4) we require initial conditions 

4>(x, t - o) = o Viefi , ( 5 ) 

and boundary conditions 

</>(£,<) Ion =0 Vt6(0,T) , 

where dtt denotes the boundary of Q, and n will denote the unit outward normal on dtl. 
The simple initial conditions (5) and homogeneous essential boundary conditions (6) are 
chosen for clarity of presentation; the method, in fact, applies to general non-homogeneous 
and non-essential boundary conditions [25], 

For the passive scalar equation (4) the convection velocity u is prescribed. The field 
u is required to be solenoidal, 

U m ,m =0 ( 7^ 
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to vanish on the domain boundary, 


( 8 ) 


and to be independent of time, 

U = u(x) . 

The requirement (8) is closely coupled to (6), and can be relaxed to include inflow and 
outflow conditions [25]. Time-dependent n will be considered in Section 4 in the context 
of the Navier-Stokes equations. 

2.2 Weak Forms 

Our numerical formulation rests on two Lagrangian variational forms, the first for 
the convection-diffusion equation (4), the second for a seemingly trivial equation which 
we shall denote the “invariance” equation. To better motivate these variational forms, 
we remark that our numerical scheme, like previous finite-element characteristic methods 
[18,19], corresponds to finite-difference methods in space-time along ^-characteristics, and 
variational methods in space. The purpose of the Lagrangian variational forms is first, 
to (Lagrangian) transform “space” into “time” as regards the convection operator, and 
second, to allow for (variational) finite or spectral element spatial discretization of the 

diffusion operator, 
a. Convection-Diffusion Equation 

To pose the variational problems we first define the classical Sobolev spaces L 2 (Q) 
and H$( ft) [26], 

L 2 (Q) — {u(x) measurable in D, J v 2 dil < oo} (10a) 

Hq(Q) = (v(x) £ I 2 (fl); v„ £ L 2 (D), q = 1,... ,d, v | 8 n = 0} : (10b) 
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In what follows we shall denote the space Hi (Cl) as Y. We shall also require the temporal 
spatial space Z defined by 


Z -{v measurable in Cl x (0 ,T), j j ( v 2 + £ (v >t ) 2 )dCldt' < oo, 

0 i 

9 — 1 

a-e. t€(0,T)} , (11) 

where T is the final time of integration. The norms associated with these spaces will be 
denoted by || . J| £3) j| . and || . j| z , respectively. 

The Lagrangian variational form of the equations (4)-(6) is then: Find 0 6 Z such 
that 0(x, t = 0) = 0, and V< e (0, T) 

— { / 00dQ} = -v J ifr'i^'idfl + J ipgdCl (12) 

S2 f i n ' 

for all 0 6 Z such that ${x y s = 0) = ^o(£.), and Vs = T — t £ (0,T) 


0>j — 0 Vx ^ fl 


(13) 


for some 0 o (x) 6 ffj(ft). Here 0(z,s) = 0(x,T - s). Equation (12) is readily derived 
from the strong form (4)-(6) by integration by parts in space and time, use of “Reynolds 
transport theorem”, and substitution of equation (13) [27]; note that the particularly 
simple form of the first term in equation (12) requires a solenoidal convecting velocity field 
u. We now suppose that the domain Cl is the image of some reference domain Cl through a 
smooth, invertible mapping A: The equation (13) for the test function 0(x,s) can 

then be written in Lagrangian form, 


/ = - 2 L( 2 L(a,s)) 

1 j;HK(a,s),s) = o 


Vs E (0,T), VaeCl 


(14a) 


/ -s = 0) = x Vq 6 H 

1 0(*, 5 = o) = 0o( £ ) \/xe ci , (146) 

where X is a standard Lagrangian spatial variable. The variable X will remain smooth 
for smooth u. 
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Note that equations (12) and ( 14) are the variational restatement of the simple phystcal 
fact that the rate of change of * in .pace-time along a characteristic of the convection 
equation is equal to the divergence of the diffusion fluxes. In our numerical scheme the 
convection initial value problem will be treated by high-order Unite-differences, and the 
diffusion boundary value problem will be treated by finite or spectral elements. 

b. Invariance Equation 

As with all Lagrangian schemes, some form of remeshing or characteristic-foot calcu- 
lation is required for the scheme to proceed for long times. A key element of our high-order 
scheme is a remesliing strategy which is both efficient and accurate; the method is based on 
an invariance procedure, in which we associate to any continuous function 5(x) £ H l 0 ( ft) 
a new function q(x,s) and related (trivial) evolution equation 


q(x,s = 0) = 5(x), g- a q(x,s) = o Vie ft • 

(15) 

Equation (15) can be written in a u-Lagrangian variational form as; Find q(x, 
that q(x, s = 0) = ?(x), and Vs € (0)^) 

s) G Z such 

- 7 - { / p ( i) + / p Um (i ‘ m ~ 0 
ds J ft J ft 

(16) 

for all p G Z such that p(x,s = 0) = po(x), and Vs G (0, S ) 


P,. “ UmP.m =0 Vxefl 

(17) 

for some po(i) € ff^(ft). Equation (17) for p(x,s) can be written in Lagrangian form, 

I ^r(a. 5 ) = — s )) Vs6(0,S), Va € ft 

l fp( 2 C(«.^)^) = ° 

(18a) 

f X{a,s = 0) = x Va G ft 

l p(x, 3 = 0) = po(x) Vx G ft , 

(185) 

with ft defined as for equation (14). 
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The equivalence of ((16), (18)) and (15) can be understood by noting that the u- 

motion of the “mesh” through (18) is cancelled in (16) by the u-convection term. Although 

it may appear overly complicated to re-introduce Lagrangian convection by the u-field 

into an equation (15) which is essentially Eulerian, this wiU allow us to accurately remesh 

variables while preserving the underlying structure of the spatial discretization. Indeed, we 

shall see that the Lagrangian simplification of (4) in (12) in exchange for the Lagrangian 

complication of (15) m (10) is equivalent to striking a better balance between At and 

Wa,,JW a 


Remark on Mow/Outdow Boundary Conditions: We note that although we have 

implicitly written the hyperbolic equation (16) as having Diriclilet boundary conditions 

on an (i.e„ q (x,s) € Z), this is, in fact, equivalent to the proper choice of no boundary 

conditions (recall u m n„ = 0), given that q(x) e //<} In treating infiow/outflow problems 

(in which we impose essential boundary conditions on * at inflow, u„A„ < 0, and natural 

boundary conditions on * at outflow, «. ». > 0), we require boundary conditions on , in 

(16) only at inflow, with this inflow value obtained through extensioii of the inflow profile 
upstream. 
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3. Discrete Equations 


3.1 Spatial Discretization 


a. Discrete Function Spaces 


We first need to generate the discrete function spaces (“nodal finite element spaces ) 
associated with the space Y = Jfj(fl), after which our method will follow directly from the 
Lagrangian variational forms of Section 2.2 and ldgli-order finite-difference discretization 
of the “time" terms in equation (12) and (16). We begin by assuming a conforming, 
non-overlapping domain decomposition [28] £>, 


K 


SI = U n 4 . 


(19) 


Jt =1 


in which the Q* are delined with respect to reference volumes 0* (e.g., a segment in JR 1 , a 
triangle or square in JR 2 , a tetrahedron or cube in JR 3 ) by an invertible elemental mapping 

(in fact, A of equation (14)), 


l k x m £ Vm 6 {1 ,..ii4 


x £ fl 


4 £ Cl k Vm £ {1 <0- 


(20a) 


(20b) 


The decomposition is geometrically conforming in that we require that the intersection 
r = dQ l n an 1 be either null, or an entire face, edge, or vertex of both Q k and fl l . 

Within each reference domain Q k ,k = {l, we Produce a polynomial space 

JP k (fl k ), a unisolvent set of N k basis points m € {l,...,d},i £ {1 >-> N )> and a set 
of associated Lagrangian interpolants Q ’ , 


Qi< k e 2P*(D*) 


Q i ' k U i ’ h ) = hi Vi,j £ N k } 2 , 


( 21 ) 
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where <5*, is the Kronecker delta symbol. For functional conformity we require that the 
lP k (Cl k ) and F‘ k be chosen such that for all k <E {l,...,A'} and for any v(Q 6 2P‘(f2‘) 
and any x 6 T _ dfl n dfl such that T is non-null there exists a w £ 2P ( (lV) such 
that w{F~ l *'(£)) = v(F~ 1>k (x)). With our nodal Lagrangian-interpolant representation 
we now make precise the particular isoparametric map F^ k of interest in our Lagrangian 
characteristic technique, 

N h 

**«*(£) = £ Vme{l d} Vk e t K} , (22) 

1=1 

in which the X ,,k are the images of the £' ,k . 

In addition to the conforming assumption on the 2P*(Q*) and £•* we require: the 
( be chosen such that those C ,k on dfl 1 are unisolvent with respect to the trace space 
of iP (ft*); the usual “local-global” coincidence condition on the X i,k (e.g., for each ( i,k 
on an edge shared by elements k and l, there must exist a unique which maps to the 
same global point), 

X_ 1,k = Xf Vi e V(j,k) e Si . (23) 

Here M is the number of global nodes in the system, and is the set of (local node, 
element) couples associated with the unique global node Xf . Note that although we 
introduce here the standard finite element global identifier <S, to describe inter-element 
connectivity and continuity, this construct serves purely for purposes of succint presenta- 
tion; in practice, our schemes are based entirely on element-local constructs (see Section 5) 
[13,14,29]. Summarizing our superscript notation: X refers to a globally defined quantity; 

X ■* to an elementally defined quantity; to a nodally/elemen tally defined quantity; 
and X' to a nodally/globally defined quantity. 

We can now define our discrete polynomial subspace Y k (X?) C #o(ft) = Fas 

n(X‘)= span ( Z <?'•*(£) V»€{1 _V} /B ) (24 a) 

(j i k ) € S , 
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where VA; E {1, . . . , A } aud Vj 6 { • • f N } 

n >.>( r \ = l $'■*(£“ 1,4 ( aO ) ( 246 ) 

g UJ lo Vx£D , 

where the set <1 -AO/- * xcludes those *' 6 ^ N) ^ which ^ ° D ^ ^ 

space n(X‘) is thus characterized by the domain decomposition V, by h % symbohc for the 
jpi (ft*), fc = i, . . . , K , and by the the global nodes £ . The basis we choose to represent 

v h G Y h is the nodal basis used to describe the space in (24) 


Vk(x) = v \ Q >,k (s) 

»=1 (> , * ) e s i 

wl = 0 * G M)',b 


(25a) 

(256) 


where the set {1, .... W/s is the complement of {1, ...,N},b in {1, and v\ v h {2L). 

For simplicity we shall write „< = o', and explicitly write ,(£<) to denote the nodal values 

of a function q not a priori in Y h (X')- 


b. Discrete Inner Products 

Having defined our spaces in Section 3.1a, we now complete the spatial discretization 
by defining the discrete inner products associated with the continuous inner products 
required in the variational forms defined in 2.2. To begin, we remark that the convection- 
diffusion equation (12) requires two inner products, the L 2 inner product and an H l 
bilinear form. Their discrete forms with respect to our space Y h are 

v * 1 ,** €(«*)’ 

i=lj=l 

£ j. ■ < 26) 


, ft g (m * ) 2 Aii ( 2 C‘ w = 

i=ii=i 
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respectively. Here and G^({) are the Jacobian and the derivative transformation 

matrix associated with the change of variables (22), 


J '*(0 = det F m , t ({) 


(28a) 


and 


a £(LHJ = 6„, , (284) 

respectively. The invariance equation (Id) requires the L * inner product (20). as well as a 
convection operator, 

WA 3 e (m *) 2 EE D** (jlW = 

i=lj=l 

X /S 

E . (29) 

* ~1 j — 1 P.f> i** 

(«>)€*, 

In (26)-(28) t refers either to exact quadrature or numerical quadrature, depending 
the particular spatial discretization chosen (see Section 5). 

Lastly, we require the linear form 


on 


.V 


W em" E r?{x!)g = 


<=i 


•- 1 (p.t)eJi ^ 

for the inhomogeneous term in (12). 
3.2 Full Discretization 


Et/.’ E £.<?'■*(£) sM£'‘(0)J •*(()</( , 


(30) 


Our fully discrete convection-difTusion equation now corresponds to insertion of the 
discrete linear and bilinear forms of (26)-(30) into the variational form (12), followed by 


discretization of the “space-time" derivative associated with the left-hand side of (12) by 
the Qth order backward-differentiation implicit finite-difference formula (6|: Find i> ' + 

= (f)>(t n+l = (n + l)At) such that 

,- l ,=0 

/?'{E A ij (X t ' n+l )<t > r ' n+1 + ^(X‘ ;n+1 )ff :n+1 > » Vi € ( 31a ) 

i= l 

^in+i _ 0 Vi € {1, • • • ■ ( 31b ) 

Here = X r> * refers to the the base (“undeformed”) discretization of ft, and thus 

the A <J (X l;n+1 ),^(*‘ !n+1 ) are known] it can also be shown from the incompressibility 
of u and (14a) that (X ,;n+1 -’) is, in fact, independent of X', and can thus be taken 
as 5*i(X ,;n+1 ). Note that any nodal variable without tilde refers to nodal values at the 
base-discretization points X' 5 *- The backward differentiation coefficients are given in [6] 
with /3, ~ 0( 1/At), and ~ 0(1). 

It remains now to find the ?*+ l - for , > 0; ^+ l - represents the value of * at 
time t n + 1 "« at the “foot” of the characteristic whose “head” at time t n+1 is at X?' [18]. 

To this end, we use our invariance procedure (16)-(18) to determine + 1 &om the 

values <^; n+1 -«, the latter being known from previous timesteps. To wit, at each timestep 

t n+l we perform the following subproblem: Vg € {1, • ■ • . Q} 

q\rn—0 __ q 

u { ” m=0 = 


Vi 6 {1,...,^} ( 32 ) 
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Vm = 0,..',qAt/As -1 ,Vi 6 


1 

As 


{* * 

(Z B*’ (2C" m+l )4> iin+1 ~ , ' m +1 - ]T B ii (2cy" m )&'> n + 1 -'-’”' 

i ~° j- o 


2 jV 

-Z Z 7 P (2£* ” p ) <^ in +1_ * :m - * 

p =0 j =zl 



B'* (X' ::m+1 )u , ! ' m + 1 


- Z B ij (x , ''" n )u y: ' m ) = 

i=o 


2 V 

- Z Z 7 P (2c* :;m ■ f )u> ” m - " 

p = 0 i —1 


(2c i; ; m+1 -2c <;:m )/As = -Z 7pM <::m - j ’ , 


>=0 


(33a) 


(336) 


(33c) 


We have introduced here an explicit third order Adams-Bashforth scheme in s (with 
coefficients 7p [6]) to solve the invariance problem. The subproblem timestep is A* < At, 
with 27 assumed integer; the notation 4> i]n refers to timestep n, while refers to 

sub-timestep m for a given t". Equations (33a) and (33c) correspond to discretization 
of (16) and (18) respectively; in conjunction with (34), (33a) and (33c) represent the 
Lagrangian requirement (14). Determination of the 4> requires application of invariance to 
the velocity as well, (33b), as the velocity enters into the convective operator D ij , (29), 
as well as the mesh update equation for £ . Note that, in practice, we do not assume 

^ (— ) = & j (2C ” ) in (33a) and (33b), as actually calculating B ij ( £ » m ) from (26) 

appears to improve stability. 

The complete scheme comprises: Lagrangian variational forms that reduce the convec- 
tion-diffusion equation to a symmetric initial value problem along characteristics; high- 
order, implicit finite-difference schemes for integration along characteristics; standard 


16 



finite-element-like discretizations for the symmetric boundary-value-problem spatial op- 
erators; mesh invariance and higli-order (explicit) time-stepping for deducing the function 
values at previous space-time points. We now turn to a summary of the characteristics of 
the method; the next section is a heuristic anticipation of the numerical experiments of 
Section 5, and of theoretical analyses treated in [25]. 

3.3 Stability, Accuracy and Computational Complexity 

a. Stability 

For the backward-differentiation scheme (31) unconditional absolute stability can be 
intuited in several ways. First, for the Q = 1 scheme, energy arguments are readily 
constructed [18]. Second, for any Q ^ 6 the region of absolute stability in the complex 
AA t plane for the model problem u t = Xu for the backward differentiation formula includes 
the entire negative real axis, implying that for diffusion partial differential equations these 
schemes are absolutely stable. As (31) is a diffusion equation along characteristics, we 
expect unconditional stability. Indeed, it is readily shown from von-Neumann analysis 
that for the simple case of Fourier spatial discretization the modulus of the backwards 
differentiation “temporal” growth factors for space-time integration of convection diffusion 
are identical to those for time integration of standard diffusion; as expected, convection 
effects only the phase of the propagator. 

As regards the application of the third order Adams Basliforth technique to the in- 
variance equations, (33), it is clear that the usual Courant condition will apply, giving a 
condition of the form As < 0{h/\u |). We note that our Lagrangian-decoupling method 
bears a resemblance to earlier subcycle proposals [20]; an essential difference is that (33) rig- 
orously and stably decouples the convection and elliptic (or Stokes) contributions, whereas 
subcycle approaches involve repeated application of the convection operator in splitting 
fashion. 
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b. Accuracy 


Convergence of <p K to <p in at a fixed time T (or in the Z norm) is achieved as 
At -+ 0 in (31), As — * 0 in (33), and as — ► Y . The error due to space-time integration 

by the backward differentiation formula in (31) is expected to be of order (A t)^. The 
ability to obtain virtually arbitrary-order accuracy constitutes a major difference between 
the current technique and past methods [18]; high-order Adams-Moulton schemes would 
not be stable, and would also require multiple Laplacian evaluations on deformed (not 
— ’ ) domains. As concerns the temporal accuracy of the Adams-Bashforth treatment of 
the invariance equation, (33), we expect 0((As) 3 ) errors. The fact that the error can 
be controlled through As and a high-order time-stepping scheme represents a significant 
improvement over earlier remeshing schemes which are typically 0{At) [19]; low-order 
remeshing methods can result in strongly accuracy-limited tiinesteps. 

The spatial errors are proportional to the error in the best-fit approximation of the 
solution 4> by the underlying discrete space Y h [14,18,25]. In particular, for fi-type finite 
element refinement [23] (A — » oo, N k fixed) we expect algebraic convergence, and for 
p-type [30,31], or spectral element refinement [14,15,32] ( K fixed, N k -» oo) we expect 
spectral convergence (e.g., exponential convergence for infinitely smooth functions). Note 
that the spatial discretization does not only effect the treatment of the diffusion operator in 
(31): spatial discretization controls the numerical diffusion and dispersion in the invariance 
procedure (33); the allowable At in (31) is (accuracy, not stability) limited by the degree 
of geometric distortion that can be represented by the spatial approximation [25,33], 

Since the characteristic procedure transforms space into time in the convection op- 
erator, a time ’-stepping scheme for (31) limited to low-order approximation would po- 
tentially preclude the use of a high-order method in the remaining spatial operators. The 

fact that the current scheme is high-order in time allows methods that are high-order in 

1 

space to once again be considered in a characteristic framework. Although it might appear 


18 



that a variational (e.g., spectral) method in space would suggest a matching variational 
(spectral) approximation in time [24], the initial-value-problem nature of the convection 
operator indicates that finite-difference approximation will be more efficient, with none 
of the attendant "complex geometry” objections typically raised as regards high-order 
finite-difference approximations in space. 

Lastly, it should be remarked that even the low-order ( Q = 1) scheme can be quite 
accurate. In particular, we note that if the invariance procedure is exact (e.g., As small), 
then the temporal error in (31) is due solely to inexact integration of the diffusion terms 
along a characteristic. This contrasts sharply with conventional upwinding schemes [21], 
in which the majority of the numerical diffusion derives from mesli-induced inaccurate 
representation of the characteristic foot. As a result, low-order forms of (31)-(34) can 
achieve high accuracy if: gradients in 4> are small, and the diffusivity u is small , gradients 
in 4> are primarily in the cross-stream direction, that is, orthogonal to u [5]; gradients in <f> 
are large in the flow direction, but the flow velocity (and hence effective spatial increment 
As =| u | At) is “small” (such as in stagnation-flow boundary layers). 


c. Computational Complexity and Generality 

The beauty of the characteristic scheme [18] is that it reduces each timestep to a 
symmetric elliptic solve (in the case of Navier-Stokes, a symmetric Stokes solve) that can 
be very rapidly solved iteratively [4,14,16,17] ; that is, the asymmetrizing and destabilizing 
effect of the convection operator is eliminated. For our particular implementation, the price 
to be paid is the solution of a "bare” (non-pressure-corrected) convection operator (33). 
The key to the computational efficiency of the scheme, as described in detail in equation 
(3) of the Introduction, is that the increase in work associated with (33), W^ tC) is more 
than compensated for by the large A t allowed in (31). In particular, we note that the work 
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for the characteristic scheme, W eh , is given by 


W ch = 


T T 

a'wz; t ' c + 


At cr 


At' 




( 35 ) 


where the first term represents (33), and the second term (31). Here A t CT is as in (3), 
a' is inflated from (3) by virtue of the additional work associated with Q > 1 in (33), 
and Af' acc is the new accuracy limit that reflects the effect of geometric distortion and the 
accuracy of the Q th order backward differentiation scheme. Although many factors must 
be considered to determine whether, at fixed accuracy, W efc is much less than W' - *', the 
potential clearly exists for great savings. In particular, the two terms in (35) are better 
balanced than in (3), and the dominant work per timestep, Waj, is effected much less 
frequently, T/At’ acc < T/At„. In summary, the stability-limiting and work intensive 
parts of the calculation have been Lagrangian-decoupled, 

A central aspect of the efficiency of the method is the invariance procedure. Pre- 
viously proposed schemes are low-order [18,19], or use a nonlinear search and solution 
algorithm to find the 4> directly. In the latter approach, [(33a), (33b)] are replaced by a 
’’nonlinear solver” to obtain the <f> directly from the X_* given by (33c) [18]. Although the 
’’nonlinear solver” may appear faster that the time-stepping approach (33), in practice the 
opposite will typically be true: first, even for low-order approximations, the ’’nonlinear 
solver” approach requires non-elemental, unstructured calculations, leading to inefficient 
vectorization and parallelization - the time-evolution solver preserves topology and struc- 
ture; second, as the ’’nonlinear solver” approach destroys discretization structure, tensor 
product representations will not be preserved, leading to orders-of-magnitude increases 
in work for high-order approximations ([14,34], see Section 5) - the time-evolution solver 
maintains an initially tensor-product representation for all ’’time” s. 

Lastly, we briefly comment on generality. Iterative solvers, unlike direct solvers, are 
typically non-robust with respect to modifications to the equations and physics, as their 
convergence (rate) is dependent not only on discrete-equation structure, but also on the 
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structure of the spectrum. It is usually the case, however, that increased physical com- 
plexity in incompressible continua enters into the “divergence of the flux tensor ’ portion 
of the equation, not the acceleration term, and thus the Lagrangian-decoupling method 
should be generally applicable to a large class of problems. 
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4. Navier-Stokes Equations 


In this section we describe the extension of the convection- diffusion scheme of Section 
3 to the full nonlinear unsteady Navier-Stokes equations (1). There are many ways in 
which to deal with the nonlinearity; we propose here a linearization technique which is 
both simple to implement and illustrative of the general procedure. 


To begin, we consider a semi-discretization in time, in which u^(x) = u m (x,n'At ) 


are known for 0 < n‘ < n, with u"' m = 0 Vn' 6 {0, 
that is, find < +1 (x), we write u m (x,t) = (x) + 6u m (x,t), 


To advance the solution, 
and choose the particular 


linearization of (1) that yields 


u m,t + p >m + + f m in fl (36 a ) 

u tl , = 0 in n (365) 

for t near t . The linearization (36) is preferrable over other possible distributions of 6u m in 
that stability is readily shown, £ f Q u m u m < 0. We now apply our Lagrangian-decoupling 
method and variational spatial operators exactly as for convection-diffusion, except that 
now the velocity appearing in the convection operator D ij is frozen at u n during the substep 
(32)-(34), and the Laplacian operator in (31) is replaced with an appropriate symmetric 
Stokes discretization [14,15]. Note that, unlike the compressible flow case, (33a) has only 
a single “characteristic” direction, n, shared by all components of the velocity. 

As regards the accuracy of the resulting scheme, it is clear that the method (36) is 
O(At) in time due to the linearization in Su m ; this does not, however, preclude use of 
a 0((At) Q ) scheme in space-time in the Lagrangian-decoupling procedure. In essence, 
the former relates to the temporal accuracy of the integration along the characteristic, 
whereas the latter relates to the spatial accuracy of the integration; note that when a 
steady-state is achieved, the 0{At) errors vanish, whereas the spatial errors remain. For 
time dependent problems, variants of (36) which are more accurate in time are readily 
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generated by appropriate liigli-order linearizations [25j. The spatial errors due to the 
Stokes discretization are standard [1,14,35,36], and will not be discussed further. 
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5. Spectral Element Implementation 


In this section we consider the spectral element method as a particular example of 
the general variational spatial discretization described in Section 3.1. We give numerical 
results which demonstrate the convergence rate of the spectral-element-based Lagrangian- 
decoupling method for the convection-diffusion of a passive scalar. We also demonstrate 
the viability of the method for the full Navier-Stokes equations, in particular as regards 
increased stability and associated computational savings. 


5.1 Description of Spectral Element Spatial Discretization 

The spectral element method is a high-order variational method for the spatial dis- 
cretization of partial differential equations [14,15,32]. As in (19), the computational domain 
is first broken up into K subdomains (spectral elements) ft*,Jk = {l,...,Jjf}. For reasons 
of efficiency (i.e., the availability of tensor product sum factorization [34]) the elements 
ft* are taken to be line segments in 1R 1 , quadrilaterals in 1R 2 , and hexahedra (bricks) in 
2R 3 , defined with respect to the reference volumes ft* =] — 1,1[ 4 in M d by an invertible 
mapping F •* (20). Although geometrically nonconforming spectral element methods have 
been developed [37,38], we restrict ourselves here to the conforming case. 

Within each reference volume ft* the discrete polynomial subspace JP k is taken to be 
the set of all polynomials of degree < N in each spatial direction [14]; the spectral element 
spatial discretization can be characterized by the discretization pair h = ( K , N). Although 
the particular choice of basis points £ 1 does not effect the error estimates, it greatly effects 
the conditioning and sparsity of the resulting set of algebraic equations, and is critical for 
the efficiency of parallel iterative solution procedures [13,29]. In order to take advantage 
of efficient sum-factorization techniques [34], the N k = (N + l) d basis points £••* for 2P* 
are taken to be the d- tensor product of the N -f 1 Gauss-Lobatto Legendre points [39]; 
the corresponding Lagrangian interpolants <?’•* of (21) are therefore the product of d N th 
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order one-dimensional Lagrangian interpolants [14]. 

In the spectral element method the discrete inner products defined in Section 3.1b 
are based on tensor-product Gauss-Lobatto Legendre integration, that is, $ is equivalent 
to Gaussian quadrature. The quadrature points are the same as the basis points; with 
N + 1 quadrature points we can integrate exactly polynomials of degree < 2N - 1. This 
particular choice of numerical integration (denoted “consistent” integration) can be shown 
to be sufficient and optimal for most problems [14,33]. In particular, consistent integra- 
tion is sufficient for complex geometries, an important consideration when considering the 
invariance equation (32), in which the base geometry may deform significantly. 

The spectral element method combines the geometric flexibilty of a low-order method 
with the rapid convergence rate associated with spectral techniques. Considering only 
spatial errors, it can be shown that the spectral element solution to the convection- diffusion 
equation (4) or to the full Navier-Stokes equations (1) converges spectrally fast to the exact 
solution for K fixed, N -> oo, with exponential convergence obtaining for locally analytic 
geometry, data and solution. This rapid convergence rate derives from the good stability 
and approximation properties of the polynomial spaces 2P‘, and the accuracy associated 
with Gauss-Lobatto Legendre quadrature and interpolation [14,39]. Moreover, the discrete 
solution suffers from minimal numerical dispersion and diffusion, a fact which is important 
in the solution of the invariance equation (16). 

Having defined our spectral element spatial discretization, we need to solve the set of 
algebraic equations (31)-(34) corresponding to the convection-diffusion of a passive scalar. 
For the implicitly treated symmetric positive-definite elliptic kernel (diffusion) in (31) we 
employ preconditioned conjugate gradient iteration or rapidly-convergent intra-element 
multigrid techniques [16,17]; for the Stokes operator of (36) we employ splitting schemes 
[11,12] or preconditioned Uzawa methods [4,14] to reduce the procedure to elliptic sub- 
problems. The key to minimal work per iteration is the use of tensor product elements, 
spaces, bases, and (consistent) quadratures, allowing for sum-factorization and efficient 
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matrix-vector product evaluations. 


We note that the solution of (31) is stable for all time steps At ( Q < 6). However, in 
order to evaluate 4> from previous timesteps at the foot of the relevant characteristics, we 
need to solve a pure convection problem for 4> and u. The details of this evaluation are given 
m (33), in which we have chosen to use an explicit third-order Adams-Bashforth integration 
scheme; the associated spectral element Courant condition is A t CT < 0{l/ K lli N 2 ) [15] 
It is important to note that our approach to solving the invariance equation (16) preserves 
the underlying tensor-product structure of the spatial discretization, and hence no costly 
re-interpolation is necessary. Indeed, if referencing were done only to the base geometry 
interpolation would require 0{I< N 2i ) operations, compared to 0(KN d + l ) for our 
tensor-product scheme [34]. 

5.2 One-dimensional results for a passive scalar 

In this section we consider the solution of the convection-diffusion equation (4) in one 
space dimension. In the first test problem we choose fl =]0, 1[, u = 0.01, u = 1, with 
<p{x,t = 0) = sm(27ra:), and periodic boundary conditions. The exact solution is given by 
- e *sin(2n(x - ut)). First, we demonstrate the temporal accuracy of the 

new Lagrangian-decoupling method. Figure 1 shows the £ 00 -error at a fixed time t = 2 as 
a function of the time step At for a fixed spatial discretization K = 2 and N = 15; this 
spatial discretization assures that the temporal error is always dominating. The results 
clearly demonstrate that the temporal error is 0((At)9), where Q is the order of the 
backward differentiation formula in (31). Since As is typically much smaller than At, 
the temporal error resulting from solving (33) will, in general, be much smaller than the 
temporal error resulting from solving (31), at least for Q < 3. We note that for pure 
convection = 0), the only temporal error is due to the integration of the invariance 
equation (32)-(34); for this smooth problem, the accuracy improves as v decreases (in 
contrast to “mesh-upwinding”). Next, we demonstrate the spatial error of the scheme. In 
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Figure 2 we plot the L°° -error at a fixed time t = 2 as a function of the polynomial degree 
JV, for fixed K - 2 elements, with Q = 3, At = 0.01, and At/ As = 100. As expected, we 
obtain exponential convergence as long as N < 10; for AT > 10 the temporal error becomes 
dominating given our particular choice of Q and At. 

As our second one-dimensional test problem we consider the steady state solution 
to (4) in ft =]0, 1[, v = 0.1, u = 1, and essential boundary conditions <p(x = 0,t) = 0, 

0 j 1 « 

<£(x = l,t) = 1. The exact solution as t -+ oo is given by <j> = 717771 . which has a 
boundary layer of thickness u at x — 1. The numerical steady solution is obtained by 
integrating the unsteady convection-diffusion equation to a time f = 10, starting from an 
initial condition <^>(x,t — 0) = x. Figure 3 shows the L°° -error as a function of the time 
step At for Q = 1,2, 3, 4, for a fixed spatial discretization K = 2 and N = 11. Note that 
in this case the time step At has to be smaller than O(u) before the expected convergence 
rate is observed. Also note that, unlike the standard Eulerian convection scheme, in which 
the steady state error is purely spatial, the error in our new scheme is contaminated by the 
temporal discretization due to the transformation of space into time. In Figure 4 we plot 
the L°° -error as a function of the polynomial degree N for fixed K = 2, Q = 3, At = 10' , 
and At/ As = 10. The results demonstrate exponential convergence, at least for N < 10; 
for N > 10 the temporal error again becomes dominating. 

We note that it is a drawback to the current method that At must be the same at all 
points on the domain, as this does not allow for local refinement near important spatial 
structures. Fortunately, the flow “through the wall” (u m n m ^ 0) associated with the 
normal boundary layer of Fig. 3 is not typical; if u m h m — 0, the characteristics scheme is, 
in some sense, self-adaptive, as Ax oc / u dt' . 

5.3 Two-dimensional results for a passive scalar 

The first two-dimensional test problem we consider is pure convection (u = 0) of 
a passive scalar <f> in Cl =]0,1[^. The boundary conditions are periodic conditions at 
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x\ = 0 and x\ — 1, and homogeneous Neumann conditions at X2 = 0 and X 2 = 1. The 
prescribed solenoidal velocity field is u i = 1, U2 = 0, and the exact solution is given by 
4>{xux 2> t) = sin{2'n{x\ — uit))cos(2nx2)- In Figure 5 we plot the error in the H 1 semi- 
norm at a fixed time f = 5 as a function of the polynomial degree N for fixed K — 4 
spectral elements, Q = 1 and At = 0.1. The results clearly demonstrates that exponential 
convergence is achieved due to the smooth nature of the solution. Note that for this pure 
convection problem (for which our new scheme is not of practical interest, as W^ ( v , is 
effectively zero), the only temporal error is the error due to integrating the invariance 
equation (32), and thus Q = 1, At = .1 is sufficient. The subtimestep As is sufficiently 
small here that the spatial error is dominating. 

As the second test problem we consider the steady state solution to the convection- 

diffusion equation (4) in the two-dimensional domain ft : xi Gj — 1/2, 1/2[,X2 6)0, 1[. The 

given velocity field is the stagnation potential-flow u\ = x\ t u 2 = — x 2 . With boundary 

conditions <f>(xi,x 2 = 0,t) = 1 and 4>(x\,x 2 — l,t) = 0 the steady state solution depends 

/•» 

only on x 2 , and is given as <f> = 1 ; . In Figure 6 we plot the error in the 

*1 ^ ® 

H semi-norm at a fixed (effectively infinite) time t = 5 as a function of the time step At 
for Q — 1,2,3 for a fixed spatial discretization consisting of K = 4 spectral elements of 
polynomial degree N = 10. Note that we impose the exact solution as essential boundary 
conditions on the whole domain boundary. The expected convergence rate 0((At )^ ) is 
clearly demonstrated. For Q = 1 we repeat the experiment, but now imposing essential 
boundary conditions only along x 2 = 0 and x 2 = 1, with homogeneous natural (Neumann) 
boundary conditions imposed along xi = -1/2 and x\ = 1/2; the results demonstrate that 
the scheme works equally well for outflow boundary conditions. In Figure 7 we plot the 
error in the semi-norm at a fixed time t = 5 as a function of the polynomial degree N for 
fixed A = 4, Q = 2 or 3, At = 10~^, and 5 < At/ As < 10. Again, the results verify the 
spectral accuracy of the scheme. 

Our stagnation flow results are different from the one-dimensional normal boundary 
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layer in two aspects. First, the boundary layer scale only as u 1 ' 2 here, as opposed to 
the atypical 1/ for normal layers. Second, and more importantly, we note that there is no 
evidence in Fig. 6 of a boundary layer threshold effect in At, as u m n m =0 and the effective 
Ax near the wall is therefore small. In most real viscous flows, in which no-slip is applied, 
U2 would in fact vary as - x \ , yielding even better accuracy than the slip U2 = -*2 case 

studied here. 

5.4 Two-dimensional Navier-Stokes results 

We now return to the unsteady Navier-Stokes system, which was the original moti- 
vation for the new scheme due to to the potential savings in computational cost. As a 
test problem we consider flow in a xi-periodic grooved channel, in which the flow is driven 
by a constant pressure gradient in the xi-direction. These first results are intended only 
to demonstrate the stability of the method, and to hint at the potential computational 

savings possible. 

The tests described below correspond to the geometry shown in Fig. 8a for a Reynolds 
number R = ^75- = 25, where / is the constant force (pressure gradient normalized by 
density), and H is the half-channel width. Only steady state results are presented. In 
all cases the spectral element mesh (X 1 -*) is based on the K — 6 elements shown in Fig. 
8a, with N = 6 in each element. In order to determine the accuracy and efficiency of 
the characteristic-decoupling scheme we first solve the problem by using a standard semi- 
implicit “Eulerian” scheme of the type (2), in which the convective term is treated with 
a third-order Adains-Bashfortli scheme. A steady-state solution (f = 4) is reached after 
834 timesteps, with the solution shown in Figures 8b and 8c. The maximum velocity, 

normalized by is found to be U mati = 1.0. 

Next, we repeat this test using the new Lagrangian-decoupling method. We use the 
first order Q = 1 scheme (only steady-state results are examined), with a timestep At 
which is 50 times larger than that for the Eulerian scheme (note that for the Eulerian 
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scheme the Navier-Stokes times tep At is chosen from stabihty considerations, At = At cr - 
for the characteristic scheme only As in (33) is limited by A t cr ). The final time t = 4 is 
now reached in 16 timesteps, giving the results shown in Figures 9a and 9b; in terms of 
computational cost, the new scheme obtained the t = 4 solution 16 times faster than the 
standard approach. As regards accuracy, the two solutions look qualitatively the same, 
however, the maximum non-dimensional velocity in the characteristic case is U maa = 0.88, 
indicating a 10% error with respect to the Eulerian scheme. Although this discrepancy 
may appear to be due to the severe distortion of the mesh, shown in Fig. 9c as a plot 
of 2L l " At/A, i tests on plane Poiseuille flow (in which only distortion effects are present) 
indicate that this degree of distortion is readily handled by the h = (K = 6, N = 6) 
approximation. Reduction of A t or increasing Q reveals that the error is, in fact temporal; 
systematic numerical convergence results for Navier-Stokes are given in [25]. 

The results of the these first Navier-Stokes calculations are by no means conclusive; 
tests to determine computational savings at fixed accuracy are currently underway [25], 
Of particular importance is understanding the dependence of At' acc and W eh on temporal 
order Q, the solution-induced deformation (in particular in boundary layers and near 
singularities), the spectral (or finite) element discretization, incomplete iteration, and the 
Stokes solver efficiency. It is clear, nevertheless, that the technique holds great promise in 
striking an efficient compromise between temporal stability and equation structure. 
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Figure Captions 

Figure 1. A plot of the discretization error || <p — 4>k ||t~ at a fixed time t = 2 as a 
function of the timestep At for a fixed (accurate) spatial discretization, when solving the 
convection-diffusion equation (4) using the new Lagrangian-decoupling method. The exact 
solution is given as <fi(x,t) = e~ v * $in(2n(x — ut) in the (periodic) domain Cl — [0,1], 
with v - 0.01 and u = 1. The results show the (dominating) temporal error of the Q = 1 
(O), <5 = 2 (□)> an(1 Q = 3 (A) backward differentiation schemes. 

Figure 2. A plot of the discretization error || — ||r,«> at a fixed time t = 2 as a function 

of the polynomial degree N in the K = 2 fixed elements, when solving the convection- 
diffusion equation (4) using the new Lagrangian-decoupling method. The timestep is set 
small to minimize temporal pollution. The exact one-dimensional periodic solution is given 
as <t>{x,t) - e - , '( 2,r ) a ‘sin( 2w(x - ut) for x € Cl - [0,1], with v - 0.01 and u = 1. Spectral 
convergence is achieved for N < 10; for N > 10 the (fixed) temporal error becomes 
dominating. 

Figure 3. A plot of the steady-state discretization error || <f> — <f>h ||r,“ as a function of the 
timestep At for a fixed (accurate) spatial discretization, when solving the convection- 
diffusion equation (4) using the new Lagrangian-decoupling method. The exact one- 
dimensional steady-state boundary layer solution is given by (f> = 4 ,/„_ ^ f°r x € Cl — [0, 1] 
(u = 0.1). The results show the (dominating) temporal error of the Q = 1 (O). Q = 2 (•). 
Q - 3 (A), and Q = 4 (A) backward differentiation schemes. 

Figure 4. A plot of the steady-state discretization error || <j> — 4 >k ||i~ as a function of 
the polynomial degree JV in K = 2 fixed elements, when solving the convection-diffusion 
equation (4) using the new Lagrangian-decoupling method. The timestep is set small 
to minimize temporal pollution. The exact one-dimensional steady-state boundary layer 
solution is given by <f> = for x G Cl = [0,1] {v = 0.1). Spectral convergence is 

achieved for N < 10; for N > 10 the (fixed) temporal error becomes dominating. 
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Figure 5. A plot of the discretization error in the H l semi-norm, | <f> — <p h |i, at a fixed 
time t = 5 as a function of the polynomial degree N in K — 4 fixed elements, when solving 
the convection equation (4) (v — 0) using the new Lagrangian-decoupling method. The 
exact two-dimensional solution is given by </>(xi,x 2 .0 = sin(2n{x\ - uit)cos(2wx2) for 
(xi,x 2 ) G 0, = [0, l] 2 , with ui = l,u 2 = 0. Periodic boundary conditions are imposed at 
xi = 0 and x\ — 1, with homogeneous Neumann conditions imposed at x 2 = 0 and x 2 = 1. 
Spectral accuracy is achieved. 


Figure 6. A plot of the steady-state discretization error in the H l semi-norm, | <f> — 
4> h |i, as a function of the timestep At for a fixed (accurate) spatial discretization, when 
solving the two-dimensional convection-diffusion equation (4) using the new Lagrangian- 
decoupling method. The exact steady-state solution is given as <fr — 1 jh for 


xi G [—1/2, 1/2], x 2 G [0,1], for an imposed stagnation flow = ^1» ^2 = ~^2- With 
imposed Dirichlet boundary conditions for <fi , the results show the (dominating) temporal 
error for backward differentiation schemes of order Q — 1 (Q), Q = 2 (□), and Q = 3 (A). 
For comparison we also plot the the error for the case Q = 1 with homogeneous Neumann 
conditions imposed along xi = ±1/2 (outflow) (•). 


Figure 7. A plot of the steady-state discretization error in the H 1 semi-norm, | 4>—4>h |l» as 

a function of the polynomial degree N in K = 4 fixed spectral elements, when solving the 

two-dimensional convection-diffusion equation (4) using the new Lagrangian-decoupling 

method. The timestep is set small to minimize temporal pollution. The exact solution is 

given by 0 = 1 — 4, ~~' 3/ay ^ for xi G [—1/2, 1/2], x 2 G [0, 1], for an imposed stagnation 
J 0 

flow ui = xi, u 2 = — x 2 . Spectral accuracy is achieved. 


Figure 8. Two-dimensional Navier-Stokes solution in a xi-periodic grooved channel at 
a Reynolds number R = = 25, in which the flow is driven by a constant pressure 

gradient in the x\ direction. Starting from rest, an approximately steady state solution 
is reached at a time t = 4. Here the numerical solution is obtained using a standard 
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(Eulerian) semi-implicit scheme described by equation (2), based on the spectral element 
mesh shown in Fig. 8a. In Fig. 8b and 8c we plot the solution in terms of streamlines and 
velocity vectors, respectively. 

Figure 9. Two-dimensional Navier-Stokes solution in a xi-periodic grooved channel at 
a Reynolds number R = = 25, in which the flow is driven by a constant pressure 

gradient in the xi direction. Starting from rest, an approximately steady state solution is 
reached at a time t — 4. Here the numerical solution is obtained using the new Lagrangian- 
decoupling method. In Fig. 9b and 9c we plot the solution in terms of streamlines and 
velocity vectors, respectively, and in Fig. 9c we indicate the mesh distortion that occurs in 
the invariance-equation integration. Note that although the characteristic method allows 
for large A t, and hence rapid convergence to the steady-state, the method does not reduce 
the parametric stiffness associated with long timescales as R — -» oo. 
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